%00000000000000000000000000000000000000000000000000000000000
%THIS FILE SIMULATES A QUASI-HYPERBOLIC WEEKLY BUFFER STOCK MODEL WITH BI-WEEKLY PAYCHECK
%00000000000000000000000000000000000000000000000000000000000

%set up environment
    clear all
%print flag (0 for no print, 1 for print)
        print_flag=1;
%set parameters
    R = 1+(0.01/52); %interest rate
    rho = 1; %risk aversion
%weekly beta distribution           
        beta=[0.90 0.93 0.97];
        n_beta=size(beta,2);
        delta=0.999;
%pay_share (percent of total income)
    pay_share = 0.7;
%temp shock distribution (spread the rest of the payshare over each week)
    %you get 0 on non-paycheck week
    %mu=average uncertain income
    %sigmat - calibrated from previous paper converting monthly to weekly
    %vol
    n_at=65;mu_at=(1-pay_share);sigmat=0.10;
    [theta,p_theta]=tauchen(n_at,mu_at,0,sigmat,3);
    p_theta=p_theta(1,:)';
%paycheck 
    n_p = 2;tp=(1:n_p); % time periods
    mu_at_pay=n_p*pay_share; %paycheck income as a fraction of total
%check that 2 week income adds up to 2
    assert(mu_at_pay + 2*mu_at==2);  
%set state space and consumption function
    max_a=5;
    n=10000; %size of grid
    a=(exp(exp(exp(linspace(0.00,log(log(log(max_a+1)+1)+1),n))-1)-1)-1)';  % set up triple exponential grid
%mp(grid points,diff temp shock values * diff paycheck values,week)
    mp = zeros(n,n_at,n_p); %all possible values of next period cash on hand (state x income x pay periods)
%declare functions    
    IUP=inline('c.^(-1/rho)','c','rho'); %inverse MU
    UP=inline('c.^(-rho)','c','rho'); %MU
    
%this is the grid of possible m_{t+1} values 
    for i = 1:n_p
        %in the 2nd week get a paycheck so in the first week next week is
        %pay week
        if i == 1
            mp(:,:,i) = a + repmat(theta',n,1)+mu_at_pay;
        else
            mp(:,:,i) = a + repmat(theta',n,1);
        end
    end
%%    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
% FIND THE STEADY STATE CONSUMPTION FUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
     %define current and next period c
        clear c m;
        c = repmat(a,1,n_p,n_beta); %we start next period c with consuming everything
        c_new = repmat(a,1,n_p,n_beta);
    %define current and next period m because we need it to extrapolate
        m=repmat(a,1,n_p,n_beta);
        m_new=repmat(a,1,n_p,n_beta);  
    
    %get rid of warning prompts 
    try
        warning('off','last')  
    catch
        warning('we cant turn the last warning off yet since it hasnt occured')
    end
    
    %deriv change       
        h=0.01;  
    %interp type
        interp_type='linear';
   
    %loop over different beta vlues     
    for b = 1:n_beta  
        diff=1000;   
         while diff > 1e-6
            %look over months to pay
            for j = 1:n_p
                %generate index (k) to use for next period consumption function
                if j==n_p
                    k=1;
                else
                    k = j+1;
                end

            %interpolate the consumption function next period conditional on
            %beta
                if m(1,k,b)==0
                    pp = interp1(m(:,k,b),c(:,k,b),interp_type,'pp');
                else
                    pp = interp1([0;m(:,k,b)],[0;c(:,k,b)],interp_type,'pp');
                end
            %calculate possible c_{t+1} values based on next period coh    
                cprime = ppval(pp,mp(:,:,j));
            %marginal utility of possible c_{t+1}
                uprime_prime = UP(cprime,rho);
                uprime_prime(uprime_prime>realmax)=realmax;
            %In the hyperbolic model the discountfactor is 
            %C'_{t+1}*beta*delta + (1-C'_{t+1})*delta
                %calculate MPC at each CoH point
                    mpc = (ppval(pp,mp(:,:,j)+h)-ppval(pp,mp(:,:,j)))./h;
                %expected marginal utility of c_{t+1}
                    e_RHS = ((mpc*beta(b)*delta +  (1-mpc)*delta).*uprime_prime)*p_theta;
            %back out what c is based on euler equation
                c_new(:,j,b) = IUP(R*e_RHS,rho);
            %calculate cash on hand
              m_new(:,j,b) = a + c_new(:,j,b);
            end
            diff = max(max(abs(c_new(:,:,b)- c(:,:,b)))); %take the largest difference between the two value functions
                  c=c_new;m=m_new;
         end
    end
    
    %make sure consumption doesn't exceed coh
        assert(max(max(max(c > m)))==0);
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%GENERATE SIMULATED DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%n=people, t=time periods
    discard=100;n=150;t=1000;T=t*n;discard=discard-1; 
    clear c_sim a_sim y_sim z_sim t_sim tax_refund 
    c_sim=zeros(t,n);a_sim=zeros(t,n);t_sim=zeros(t,n);y_sim=zeros(t,n);z_sim=zeros(t,n);

    %generate random sequences
        cutoff=cumsum(p_theta',2); %calculate the discrete cutoff points in the distribution of income
        cutoff_beta=[1/n_beta:1/n_beta:1]; %calculate the discrete cutoff points in the distribution of beta
        random_y=zeros(T,1);
        rand('seed',1234);
        epsilon=rand(T,1);
        epsilon_beta=rand(n,1);

     %assigns the random components based on the discretized distribution
      %the idea is to take the N(0,1) draw and see where it falls in the cumulative
      %dist. Then assign a discrete value to it
        idx1=(n_at+1)-sum(repmat(cutoff,T,1) >= repmat(epsilon,1,n_at),2); %look how far we are in the CDF       
        idx_beta=(n_beta+1)-sum(repmat(cutoff_beta,n,1) >= repmat(epsilon_beta,1,n_beta),2); %look how far we are in the CDF
            idx_beta(idx_beta>n_beta)=n_beta;
        random_y=reshape(theta(idx1),t,n); %reshape the random shocks back to matrix form
               
        %create simulated y
            y_sim =  random_y + repmat(ismember(mod(1:t,n_p),0)',1,n).*mu_at_pay; %create the simulated income values

    %start with asset value of 1
            z_sim(1,:)=1;t_sim(1,:)=1;a_sim(1,:)=0;

%generate history of CoH and consumption            
%we have a different consumption function for each period
for b = 1:n_beta
    b;
    beta_index = (idx_beta == b);
        for i=2:t
               period_index=mean(t_sim(i-1,beta_index));
            
                %Step 1: What was consumption last period? This is based on the
                %starting value of CoH in period 1
                    pp = interp1([0;m(:,period_index,b)],[0;c(:,period_index,b)],interp_type,'pp');  
                    c_sim(i-1,beta_index)=ppval(pp,z_sim(i-1,beta_index));

                %Step 2: Calculate CoH
                    a_sim(i-1,beta_index) = z_sim(i-1,beta_index) - c_sim(i-1,beta_index);
                    z_sim(i,beta_index) =  y_sim(i,beta_index)  + a_sim(i-1,beta_index);    

                %Step 3: Iterate time period 
                    if t_sim(i-1,beta_index)==n_p
                        t_sim(i,beta_index)=1;
                    else
                        t_sim(i,beta_index)=t_sim(i-1,beta_index)+ 1;
                    end
        end
end

    c_sim(end,:)=[];z_sim(end,:)=[];y_sim(end,:)=[];a_sim(end,:)=[];t_sim(end,:)=[];t=t-1;
    c_sim(1:discard,:)=[];z_sim(1:discard,:)=[];y_sim(1:discard,:)=[];a_sim(1:discard,:)=[];t_sim(1:discard,:)=[];t=t-discard;

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%
%Calculate log deviations from average [PAY PERIOD]
%%%%%%%%%%%%%%%%%%%%%%%%%%
clear log_avg_liq rel_liq cg cg_pay xx2 ind yy2 cg cg_pay reshape_cg reshape_term reshape_liq
    %calcualte average liquidity conditional on pay or non-pay period
        pay_period_index=(t_sim==2); %index for pay period
        no_payperiod_index=(t_sim==1);
        
        log_avg_liq_nopay = log(mean(a_sim(no_payperiod_index(:,1),:))); %avg liquidity for each person in pervious period
            prev_coh = [zeros(1,n);a_sim(1:end-1,:)];
                %make 0 very small numbers
                prev_coh(prev_coh<=0)=eps;
            rel_liq = log(prev_coh(pay_period_index(:,1),:)) - repmat(log_avg_liq_nopay,size(a_sim(period_index(:,1),:),1),1); %rel liq for each person
                                           
    %calcualte consumption growth
        cg =log(c_sim(2:end,:))-log(c_sim(1:end-1,:)); %consumption grown from t to t+1
        cg_pay=cg(pay_period_index(1:end,1),:); %capture all the pay_periods 

    %set colors and line types
         C = {[0 0 0.5],[0 0.5 0],[0.5 0.5 0.5]}; % Cell array of colors.
         type = {'-','--','-.'}; % Cell array of colors.
        
    clf
    hold on
    for b = [1 2 3]
        beta_index = (idx_beta == b);
        reshape_term = size(rel_liq(:,beta_index),1) * size(rel_liq(:,beta_index),2);
        reshape_liq = reshape(rel_liq(:,beta_index),reshape_term,1);
        reshape_cg = reshape(cg_pay(:,beta_index),reshape_term,1);  
        yy2 = smooth(reshape_liq,reshape_cg,0.90,'lowess');    
        [xx2,ind] = sort(reshape_liq);
        plot(xx2,yy2(ind),'color','k','LineStyle',type{b},'LineWidth',2);
    end  
    
    %PLOT FIGURE
        hold off
        axis([-2,2,-0.21,0.11]);
        ylabel('Expected consumption growth');
        xlabel('Relative liquidity');
        legendInfo=sprintf('\\beta=%.2f',beta(1));
        legendInfo2=sprintf('\\beta=%.2f',beta(2));
        legendInfo3=sprintf('\\beta=%.2f',beta(3));
        h_legend=legend(legendInfo,legendInfo2,legendInfo3,'Location','southoutside','Orientation','horizontal');
        h_legend.FontSize=12;
        if print_flag==1 
            print -depsc 'fig\qh_gc_plot_pay.eps'
        end    
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%
%Calculate log deviations from average [NON-PAY PERIOD]
%%%%%%%%%%%%%%%%%%%%%%%%%%
clear log_avg_liq rel_liq cg cg_pay xx2 ind yy2 cg cg_pay reshape_cg reshape_term reshape_liq
    %calcualte average liquidity conditional on pay or non-pay period
        log_avg_liq_pay = log(mean(a_sim(pay_period_index(:,1),:))); %avg liquidity for each person in pervious period
            prev_coh = [zeros(1,n);a_sim(1:end-1,:)];
                 prev_coh(prev_coh<=0)=eps;
            rel_liq = log(prev_coh(no_payperiod_index(:,1),:)) - repmat(log_avg_liq_pay,size(a_sim(period_index(:,1),:),1),1); %rel liq for each person

    %calcualte consumption growth
        cg = log(c_sim(2:end,:))-log(c_sim(1:end-1,:)); %consumption grown from t to t+1
        %add in a row of zeros so the dimensions are correct
            cg=[cg;zeros(1,n)];
        cg_pay=cg(no_payperiod_index(1:end,1),:); %capture all the non_pay_periods 
        
clf
hold on
for b = [1 2 3]
    beta_index = (idx_beta == b);
    reshape_term = size(rel_liq(:,beta_index),1) * size(rel_liq(:,beta_index),2);
    reshape_liq = reshape(rel_liq(:,beta_index),reshape_term,1);
    reshape_cg = reshape(cg_pay(:,beta_index),reshape_term,1);  
    yy2 = smooth(reshape_liq,reshape_cg,0.90,'lowess');    
    [xx2,ind] = sort(reshape_liq);
    plot(xx2,yy2(ind),'color','k','LineStyle',type{b},'LineWidth',2);
end    

    axis([-0.61,0.61,-0.2,0.4]) 
    ylabel('Expected consumption growth')
    xlabel('Relative liquidity')
    legendInfo=sprintf('\\beta=%.2f',beta(1));
    legendInfo2=sprintf('\\beta=%.2f',beta(2));
    legendInfo3=sprintf('\\beta=%.2f',beta(3));
    h_legend=legend(legendInfo,legendInfo2,legendInfo3,'Location','southoutside','Orientation','horizontal');
    h_legend.FontSize=12;
    if print_flag==1 
        print -depsc 'fig\qh_gc_plot_nopay.eps'        
    end